Method for improving small disturbance stability after double-fed unit gets access to the system

ABSTRACT

A method for improving system small disturbance stability after double-fed unit gets access to the system belongs to the field of electric power system operation and control technology. A sensitivity analysis is adopted to optimize parameter, through making sensitivity analysis on the non-ideal dominant mode happens to the system to find out several nonzero elements that are most sensitive to this mode in system matrix; elements of state matrix is adopted to replace the elements of system matrix to make analysis so as to find out the most relevant parameter set; setting parameters change in the interval to observe track for the change of eigenvalues of corresponding mode and then balancing and optimizing system parameters comprehensively according to the change of eigenvalues. Without adding other control means, the present invention can improve dominant modal damping caused by selecting improper controller parameters or system parameters after double-fed unit gets access to the system without increasing cost; as this method is also highly targeted, exhaustive efforts for all the adjustable parameters of the system can be avoided, which not only greatly decreases workload, but also improves computational efficiency, so that it is very instructive.

This is a U.S. national stage application of PCT Application No. PCT/CN2014/000347 under 35 U.S.C. 371, filed Mar. 28, 2014 in Chinese, claiming priority benefit of Chinese Application No. 201310109236.4, filed Mar. 29, 2013, which is hereby incorporated by reference.

TECHNICAL FIELD

The present invention relates to a method for improving small disturbance stability after double-fed unit gets access to a system, especially a method with sensitivity analysis to improve disturbance stability after double-fed unit gets access to a system, which belongs to the field of operation and control technology of an electric power system.

BACKGROUND TECHNOLOGY

With reserves of fossil fuel decrease, environmental pollution and energy shortage have been predominant and highly concerned. In order to solve the above problem, China is now striving to develop new energies such as wind power, and wind power generation has stepped into steady development after rapid growth for many years. Although the problem of wind power grid connection has been eased to a certain extent in recent years, after wind power gets access into the system, it will still cause a series of problems, which limits the further development of the wind power. Among the problems, small disturbance stability problem generated after grid connection of wind power is especially acute. At present, the mainstream model of wind turbines in China is double-fed induction generator. After getting access to the system, the generators will lead to weak damping mode, reduce system small disturbance stability margin and increase the risks for the system to operate normally. Therefore, it is very urgent to research how to improve small disturbance stability after double-fed unit gets access to the system.

Traditional methods for improving small disturbance stability after wind turbines get access to the system mainly include adopting some algorithm to optimize overall parameters or adding PSS device. The former method is of a large workload with modest effect and it cannot effectively control certain several dominant modes; the latter one is of obvious effect, but the parameter adjustment is complex and the investment cost will also rise. For the method of using sensitivity analysis to optimize system parameters mentioned in the present invention, sensitivity analysis can be made on the weak damping mode and relevant parameters can be comprehensively selected through analyzing change of eigenvalues according to analysis results of system characteristic matrixes. In this way, for a series of unstable modes or weak damping modes generated for selecting parameters unreasonably, relevant modal damping can be improved without adding other devices. Besides, this method also can greatly decrease the workload generated by unnecessary try, so that it is instructive for how to set system parameters.

SUMMARY OF THE INVENTION

The present invention, first, obtains the system matrix A′ by using eigenvalue analysis after double-fed unit gets access to the system, then through eigenvalue analysis finds out the dominant mode that requires special attention (the dominant mode includes unstable modes or weak damping modes), determines sensitivity of the elements in matrix A′ of the above modes and selects one to two nonzero elements with the highest sensitivity; then sequentially changes the changeable parameters of the elements with the best sensitivity (the specific elements need to be determined according to the sensitivity analysis results) within a certain interval (the interval is related to the elements that need to be changed; if the elements are electric parameters, the interval is the parameter value range of the actual electric components; if the elements are control parameters, the interval is the upper and lower limits of the specific controller parameters; if no upper and lower limits are provided for the parameters, the interval can be chosen freely under the condition that the lower limit is not negative) and observes the trend for eigenvalues of the dominant mode to determine optimal value of all parameters.

Technical solution of the present invention refers to a method for improving small disturbance stability after double-fed unit gets access to the system through adopting sensitivity analysis to optimize controller and system parameters, comprising the following steps:

-   Step 1: building complete mathematical models for double-fed unit,     which mainly include aerodynamic model, generator model, mechanical     model and control system model, listing a system state equation and     an output equation and then building small disturbance mathematical     model Δ{dot over (x)}=A′Δx by integrating system power flow equation     after double-fed unit gets access to the system. The system power     flow equation is current technology, where the specific format     relates to the selection of the system output variable and the     object is to constitute simultaneous equations with the output     equation in order to offset the output variable and set up the     relationship between the state variable and the input variable. -   Step 2: calculating a left modal matrix ψ and a right modal matrix φ     of a matrix A′, determining the sensitivity of unstable modes or     weak damping modes to the matrix A′ with the formula

$\frac{\partial\lambda_{i}}{\partial a_{k\; j}} = {\psi_{ik}\phi_{ji}}$ and finding out one to two nonzero elements a′_(ij) with the highest sensitivity in the matrix; analysis indicates that at the low and middle frequency band concerned by small disturbance stability of the system, the difference between eigenvalue of the matrix A′ and that of a state matrix A is not very large and as expression of A′ is very complex, element a_(ij) in A is used to make sensitivity analysis instead of the element a′_(ij) in A′. The formula is the sensitivity calculation formula where ∂ is a partial derivative, λi is the i-th mode, a′_(kj), ψ_(kj), φ_(kj) are elements at the k-th row and j-th column of the system matrix, the left mode matrix and right mode matrix, respectively.

-   Step 3: for controller or system parameters that can be set in     a_(ij), changing value of these parameters within a certain     interval, observing steady-state value of the variable required     while calculating eigenvalue of matrix A′ in simulation results,     then putting the steady-state value in A′ to solve the eigenvalues     that correspond to each group of parameters and drawing a chart on     the change track of eigenvalues. If eigenvalues are disperse, a part     of the overlapped eigenvalues need to be locally enlarged to observe     the trend for eigenvalues of dominant modes of the system. -   Step 4: if there are other parameters that can be set in a_(ij),     then repeating Step 3. -   Step 5: comprehensively analyzing the chart on modal eigenvalues     change with the track of parameters in Step 4, adjusting the     parameters in Step 4, then selecting appropriate parameter     combination upon comparison, with which both dominant modal damping     and small disturbance stability margin of the system can be     obviously improved.

The system matrix A′ illustrated in the Step 1 is built in the following methods:

Selecting appropriate state variable, input variable and output variable; the state equation, output equation and power flow equation of the system can be expressed in the following general forms: {dot over (x)}=f(x,u) y=g(x,u) y=h(x,u)  (1)

Linearizing Equation (1) at steady-state operating point, it can be concluded as: Δ{dot over (x)}=AΔx+BΔu Δy=CΔx+DΔu Δy=EΔx+FΔu  (2) Wherein,

$\begin{matrix} {A = {{\begin{bmatrix} \frac{\partial f_{1}}{\partial x_{1}} & \ldots & \frac{\partial f_{1}}{\partial x_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial f_{n}}{\partial x_{1}} & \ldots & \frac{\partial f_{n}}{\partial x_{n}} \end{bmatrix} B} = {{\begin{bmatrix} \frac{\partial f_{1}}{\partial u_{1}} & \ldots & \frac{\partial f_{1}}{\partial u_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial f_{n}}{\partial u_{1}} & \ldots & \frac{\partial f_{n}}{\partial u_{n}} \end{bmatrix} C} = {\quad{{\left\lbrack \begin{matrix} \frac{\partial g_{1}}{\partial x_{1}} & \ldots & \frac{\partial g_{1}}{\partial x_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial g_{n}}{\partial x_{1}} & \ldots & \frac{\partial g_{n}}{\partial x_{n}} \end{matrix} \right\rbrack D} = {{\left\lbrack \begin{matrix} \frac{\partial g_{1}}{\partial u_{1}} & \ldots & \frac{\partial g_{1}}{\partial u_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial g_{n}}{\partial u_{1}} & \ldots & \frac{\partial g_{n}}{\partial u_{n}} \end{matrix} \right\rbrack E} = {\quad{{\left\lbrack \begin{matrix} \frac{\partial h_{1}}{\partial x_{1}} & \ldots & \frac{\partial h_{1}}{\partial x_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial h_{n}}{\partial x_{1}} & \ldots & \frac{\partial h_{n}}{\partial x_{n}} \end{matrix} \right\rbrack F} = \left\lbrack \begin{matrix} \frac{\partial h_{1}}{\partial u_{1}} & \ldots & \frac{\partial h_{1}}{\partial u_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial h_{n}}{\partial u_{1}} & \ldots & \frac{\partial h_{n}}{\partial u_{n}} \end{matrix} \right\rbrack}}}}}}}} & (3) \end{matrix}$

Joining with Equation (2), it can be concluded as: Δ{dot over (x)}=A′Δx  (4) wherein, A′=A+B(F−D)⁻¹(C−E)  (5)

In Step 2, on the basis of the system matrix A′ obtained from Step 1, finding out the mode λ_(i),i=1,2, . . . ,m that needs to be focused, wherein, “m” refers to the number of unstable or weak damping modes. Then finding out left modal matrix ψ′ and right modal matrix φ′ of matrix A′:

For left modal matrix ψ′: ψ′=[ψ′₁ ^(T)ψ′₂ ^(T) . . . ψ′_(n) ^(T)]^(T)  (6) wherein, ψ′_(i)A′=λ_(i)ψ′_(i),i=1,2, . . . ,n  (7) “n” is the number of state variables;

For right modal matrix φ′: φ′=[φ′₁φ′₂ . . . φ′_(n)]  (8)

-   -   wherein,         A′φ′_(i)=λ_(i)φ′_(i),i=1,2, . . . ,n  (9)

The sensitivity of eigenvalue λ_(i) to the element of A′ can be expressed as:

$\begin{matrix} {\frac{\partial\lambda_{i}}{\partial a_{kj}^{\prime}} = {{\psi_{i}^{\prime}\frac{\partial A^{\prime}}{\partial a_{kj}^{\prime}}\phi_{i}^{\prime}} = {\psi_{ik}^{\prime}\phi_{ji}^{\prime}}}} & (10) \end{matrix}$

The sensitivity of eigenvalue λ_(i) to a′_(kj) quantizes the change scope of λ_(i) when a′_(kj) changes, namely, when a′_(kj) changes,

$\frac{\partial\lambda_{i}}{\partial a_{kj}^{\prime}}$ is larger, λ_(i) changes more obviously. Thus, after obtaining eigenvalues of A′, for the unstable mode, weak damping evanescent mode and weak damping ratio oscillation mode that may influence small disturbance stability of the system directly, the nonzero element that is most sensitive to this mode can be found according to the above method. However, it can be concluded from Equation (5) that A′ is obtained through a series of operations like matrix multiplication and inversion, but specific expression is hard to get. Equation (5) can be adapted to: A′=A+A _(other)  (11) wherein, A _(other) =B(F−D)⁻¹(C−E)  (12)

From Equation (11), it can be visually seen that the state matrix A is a component of the system matrix A′, so that the corresponding elements of A also exist in A′. Then after obtaining

$\frac{\partial\lambda_{i}}{\partial a_{kj}^{\prime}},$ i=1,2, . . . ,m and finding out the highly sensitive set of the elements {a′_(kj)} that are correspond to mode λ_(i),i=1,2, . . . ,m, analyzing with the component {a_(kj)} of the state matrix that shares the same code with the elements in {a′_(kj)}.

The reason for finding nonzero element is that: for a system with a fixed structure, the structure of its system matrix A′ is also fixed. If a_(kj)=0, no matter how to change parameters, a_(kj) remains unchanged.

In the Step 3, on the basis of finding out the set {a_(kj)} of elements with high sensitivity in Step 2, finding out the adjustable controller parameters or system parameters in {a_(kj)}. Wherein, k_(i),i=1,2, . . . ,t, “t” refers to the number of adjustable variables in {a_(kj)}. Setting k₁ as an example, letting k₁ change within the set interval [a,b], selecting several parameter nodes within this interval and then observing steady-state value of the variables required while calculating eigenvalues of A′. If the steady-state value basically remains unchanged, then selecting small step size Δk_(s) and cycle calculating eigenvalues of A′ while k_(i),i=1,2, . . . ,t is changing; if steady-state value changes obviously, then selecting large step size Δk_(i), taking down each simulation steady-state value and calculating eigenvalue of A′.

After working out the eigenvalues of each group of A′ while k₁ is changing, arranging these eigenvalues in a certain order (such as λ₁,λ₂, . . . ,λ_(n)) and drawing the chart on the track for changes of eigenvalues. For similar eigenvalues, please use different marks (“*”, “Δ”, “⋆”, “∘”) while drawing to avoid confusion. Then selecting optimal {circumflex over (k)}₁ from the chart on the track for changes. What needs to be noted is that k₁ might be related to certain several modes, so that it needs to balance changes of other dominant modes while selecting {circumflex over (k)}₁.

In Step 4, on the basis of selecting {circumflex over (k)}₁, repeating Step 3 for other adjustable parameters k_(i),i=2, . . . ,t, until all the adjustable parameters are set.

In Step 5, comprehensively compare analysis results of the eigenvalues of optimal parameter set {circumflex over (k)}_(i),i=1,2, . . . ,t and original parameter set k_(i),i=1,2, . . . ,t. Upon analysis, modal damping of the system can be greatly improved and small disturbance stability margin of the system can be intensified after optimizing parameters with sensitivity analysis.

The present invention is to improve the dominant modal damping generated by selecting improper controller parameters or system parameters after double-fed unit gets access to the system without adding other control means. Traditional methods for optimizing parameters are not highly targeted and for a large-scale system, traditional methods are of large workload, low efficiency and relatively blindness. For the method mentioned in the present invention, sensitivity analysis can be made on the non-ideal dominant mode λ_(i),i=1,2, . . . ,m happened to this system directly to find out the set of state matrix elements with highest sensitivity {a_(kj)} of this mode, so as to find out the set of most relevant parameters k_(i),i=1,2, . . . ,t. Then, it only needs to optimize the parameters in this set, which will greatly reduce workload. Parameter k_(i),i=1,2, . . . ,t can be analyzed through optimizing one by one, through which the track for the change of eigenvalue that corresponds to mode λ_(i),i=1,2, . . . ,m while each parameter is changing can be obtained visually. Finally, comprehensively analyzing change tracks of all the parameters can confirm a group of appropriate parameter set {circumflex over (k)}_(i),i=1,2, . . . ,t. This group of parameters can be used to improve damping characteristic of this system obviously. This method has not increased cost, avoided BruteForee/exhaustive efforts on all adjustable parameters of the system and improved working efficiency, which is instructive for how to set system parameters from improving small disturbance stability.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a system structure drawing for double-fed wind turbine-infinite system under the PSCAD/EMTDC of the present invention;

FIGS. 2a and 2b are the tendency charts for eigenvalues of the system mode while T_(issc) is changing of the present invention;

FIGS. 3a and 3b are the tendency charts for eigenvalues of the system mode while K_(pssc) is changing of the present invention;

FIGS. 4a and 4b are the tendency chart for eigenvalues of the system mode while L is changing of the present invention; and

FIG. 5 is a flow chart of the method to improve disturbance stability after double-fed unit gets access to a system according to the present invention.

SPECIFIC EMBODIMENTS OF THE PRESENT INVENTION

With reference to the drawings and embodiments, the present invention is illustrated in details below. It should be noted that the following illustrations are only examples, but not to limit the range and applications of the present invention.

The present invention adopts the model that is commonly used in researching characteristics of grid-connected operation of double-fed unit in PSCAD/EMTDC: DFIG_2010_11 (WIND FARM Vector Controlled Doubly-Fed Induction Generator) as a template to make the following improvements on this model:

-   1) adjusting reference frequency of the system from 60 Hz to 50 Hz; -   2) removing “crowbar” circuit to make it suitable for small     disturbance stability analysis after double-fed unit connects to the     grid.

This model is double-fed wind turbine—infinite system. Correctness and practicability of the present invention will be verified with this model below. The system structure of this model shall refer to FIG. 1.

In Step 1, system matrix A′ is formed in the following process:

First, building a state equation and an output equation of the double-fed unit. The system state equation is expressed as:

$\begin{matrix} {\mspace{79mu}{{{\overset{.}{\omega}}_{m}^{\prime} = {\frac{1}{T_{j}}\left( {\frac{C_{p}A_{blade}\rho\; v^{3}}{2S_{B}\omega_{m}^{\prime}} - \frac{\gamma\; u_{s}L_{m}i_{qr}}{\alpha\; L_{s}S_{B}} - {D\;\omega_{m}^{\prime}}} \right)}}\mspace{79mu}{{\overset{.}{\omega}}_{ref}^{\prime} = {{{- \frac{1}{T_{v}}}{\overset{.}{\omega}}_{ref}^{\prime}} + {\frac{1}{\lambda_{eq}\; T_{v}}v}}}{{\overset{.}{i}}_{qr} = {{\left( {T_{i\; 2} - \frac{{DK}_{p\; 2}}{T_{j}}} \right)\omega_{m}} + \frac{K_{p\; 2}{kv}^{3}}{T_{j}\omega_{m}} - {\left( {T_{i\; 2} - \frac{K_{p\; 2}}{T_{v}}} \right)\omega_{ref}} - {\frac{K_{p\; 2}\gamma\; u_{s0}L_{m}}{\alpha\; T_{j}S_{B}L_{s}}i_{qr}} - {\frac{K_{p\; 2}}{\lambda_{eq}T_{v}}v}}}\mspace{79mu}{{\overset{.}{i}}_{dr} = {\frac{{\alpha T}_{il}L_{s}}{{\alpha\; L_{s}} + {K_{pl}\;\gamma\; u_{s}L_{m}}}\left\lbrack {Q_{ref} - {\frac{\gamma\; u_{s}}{\; L_{s}}\left( {\frac{i_{dr}L_{m}}{\alpha} - \frac{u_{s}}{\omega_{s}}} \right)}} \right\rbrack}}\mspace{79mu}{{\overset{.}{x}}_{1} = {\frac{1}{T_{x}}\left( {u_{dc} - x_{1}} \right)}}\mspace{79mu}{{\overset{.}{i}}_{1{dref}} = {{\left( {\frac{K_{pcvc}}{T_{x}} - T_{icvc}} \right)x_{1}} - {\frac{K_{pcvc}}{T_{x}}u_{dc}} + {T_{icvc}u_{dcref}}}}{{\overset{.}{u}}_{1{dref}} = {{{K_{pssc}\left( {\frac{K_{pcvc}}{T_{x}} - T_{icvc}} \right)}x_{1}} + {T_{issc}i_{1{dref}}} - {\frac{K_{pssc}}{L}u_{1{dref}}} + {\left( {\frac{K_{pssc}R}{L} - T_{issc}} \right)i_{1d}} - {\frac{K_{pssc}K_{pcvc}}{T_{x}}u_{dc}} + {K_{pssc}T_{icvc}u_{dcref}}}}\mspace{79mu}{{\overset{.}{i}}_{1d} = {{\frac{1}{L}u_{1{dref}}} - {\frac{R}{L}i_{1d}}}}\mspace{79mu}{{\overset{.}{u}}_{1{qref}} = {{\left( {\frac{K_{pssc}R}{L} - T_{issc}} \right)i_{1q}} - {\frac{K_{pssc}}{L}u_{1{qref}}} + {T_{issc}i_{1{qref}}}}}\mspace{79mu}{{\overset{.}{i}}_{1q} = {{\frac{1}{L}u_{1{qref}}} - {\frac{R}{L}i_{1q}}}}\mspace{79mu}{{\overset{.}{u}}_{dc} = \frac{\gamma\;{u_{s}\left\lbrack {i_{1\; d} - {\left( {1 - \omega_{m}} \right)\frac{K_{m}}{\alpha\; L_{s}}i_{qr}}} \right\rbrack}}{{Cu}_{dc}}}}} & (13) \end{matrix}$

In Equation (13), T_(j) is an inertia time constant of the generator; ω_(s) and ω′_(m) refer to synchronous speed and mechanical speed of the generator; C_(p) is a rotor power coefficient; A_(blade) is a blade area; ρ is the air density; ν is wind speed; S_(B) is benchmark capacity of the system; γ is coordinate transformation coefficient; α is stator-rotor turns ratio; u_(s) is phase voltage amplitude of generator stator; L_(s) and L_(m) refer to generator stator self-inductance, stator-rotor mutual inductance and generator rotator self-inductance; D is generator damping coefficient; ω′_(ref) is reference wind speed used to realize MPT; λ_(eq) is equivalent tip speed ratio; T_(v) and T_(x) refer to time constants during inertia loop; i_(dr) and i_(qr) refer to the current of rotator at axis “d” and axis “q”; i_(1d) and i_(1q) refer to the current of grid-side converter at axis “d” and axis “q”; R and L refer to converter resistance and reactance; K_(p1), T_(i1), K_(p2), T_(i2), K_(pcvc), T_(icvc), K_(pssc) and T_(issc) are PI controller parameters; u_(dcref), u_(dc) and x₁ refer to reference value, actual value and measured value of DC voltage; u_(1dref) and u_(1qref) refer to inner ring PI output of grid-side controller; i_(1dref) and i_(1qref) refer to outer ring PI output of grid-side controller; C is a DC capacitor.

The output equation is expressed as:

$\begin{matrix} {{P_{g} = {1.5{u_{s}\left( {\frac{L_{m}i_{qr}}{\alpha\; L_{s}} - i_{1d}} \right)}}}Q_{g} = {{- 1.5}{u_{s}\left( {\frac{{au}_{s} - {\omega_{s}L_{m}i_{dr}}}{{\alpha\omega}_{s}L_{s}} - i_{1q}} \right)}}} & (14) \end{matrix}$

The power flow equation of the system is expressed as: P _(g) =U ² G _(line) −UE(G _(line) cos θ+B _(line) sin θ) Q _(g) =−U ² B _(line) −UE(G _(line) sin θ−B _(line) cos θ)  (15)

In Equation (15), U and θ refer to an effective value and a phase angle of a high-side voltage of a transformer; E is an effective value of a line voltage of an infinite electric power bus; G_(line) and B_(line) refer to overall equivalent conductance and susceptance of the transformer and line.

Referring to the method in (2)-(3), linearizing and organizing Equations (13)-(15) at steady-state operation point, it can be concluded as: A′=A+B(F−D)⁻¹(C−E)  (16) wherein, state variable Δx, input variable Δu and output variable Δy can be expressed as: Δx=[Δω*_(m)Δω*_(ref)Δi_(dr)Δi_(qr)Δx₁Δi_(1dref)Δu_(1dref)Δi_(1d)Δu_(1dref)Δi_(1q)Δu_(dc)]^(T) Δu=[Δu_(s)Δθ]^(T) Δy=[ΔP_(g)ΔQ_(g)]^(T)  (17)

The matrix that corresponds to Equation (16) is shown as:

$A = \overset{\begin{matrix} {{\mspace{101mu}\mspace{31mu}\;}{\Delta\omega}_{m}^{*}} & {\Delta\omega}_{ref}^{*} & {\;{\Delta i}_{dr}\mspace{14mu}} & {\Delta i}_{qr} & {\Delta x}_{1\mspace{11mu}} & {\Delta i}_{1{dref}} & {\Delta u}_{1{dref}} & {\Delta i}_{1d} & {\Delta u}_{1{dref}} & {\Delta i}_{1q} & {\Delta u}_{dc} & \; \end{matrix}}{\begin{bmatrix} {\Delta{\overset{.}{\omega}}_{m}^{*}} & a_{0101} & 0 & 0 & a_{0104} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {\Delta{\overset{.}{\omega}}_{ref}^{*}} & 0 & a_{0202} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {\Delta{\overset{.}{i}}_{dr}} & 0 & 0 & a_{0303} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {\Delta{\overset{.}{i}}_{qr}} & a_{0401} & a_{0402} & 0 & a_{0404} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {\Delta{\overset{.}{x}}_{1}} & 0 & 0 & 0 & 0 & a_{0505} & 0 & 0 & 0 & 0 & 0 & a_{0511} \\ {\Delta{\overset{.}{i}}_{1{dref}}} & 0 & 0 & 0 & 0 & a_{0605} & 0 & 0 & 0 & 0 & 0 & a_{0611} \\ {\Delta{\overset{.}{u}}_{1{dref}}} & 0 & 0 & 0 & 0 & a_{0705} & a_{0706} & a_{0707} & a_{0708} & 0 & 0 & a_{0711} \\ {\Delta{\overset{.}{i}}_{1d}} & 0 & 0 & 0 & 0 & 0 & 0 & a_{0807} & a_{0808} & 0 & 0 & 0 \\ {\Delta{\overset{.}{u}}_{1{dref}}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & a_{0909} & a_{0910} & 0 \\ {\Delta{\overset{.}{i}}_{1q}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & a_{1009} & a_{1010} & 0 \\ {\Delta{\overset{.}{u}}_{dc}} & a_{1101} & 0 & 0 & a_{1104} & 0 & 0 & 0 & a_{1108} & 0 & 0 & a_{1111} \end{bmatrix}}$

The subscript “₀” represents steady-state value of relevant variable and superscript “*” represents per unit value. The corresponding matrix elements are shown as:

${a_{0101} = {{- \frac{1}{T_{j}}}\left( {\frac{{kv}^{3}}{\omega_{m\; 0}^{2}} + D} \right)}},{a_{0104} = {- \frac{\gamma\; u_{s\; 0}L_{m}}{\alpha\; T_{j}S_{B}L_{s}}}},{a_{0202} = {- \frac{1}{T_{v}}}},{a_{0303} = {- \frac{\gamma\; T_{i\; 1}u_{s\; 0}L_{m}}{{\alpha\; L_{s}} + {K_{p\; 1}\gamma\; u_{s\; 0}L_{m}}}}},{a_{0401} = {T_{i\; 2} - \frac{K_{p\; 2}D}{T_{j}} - \frac{K_{p\; 2}{kv}^{3}}{T_{j}\omega_{m\; 0}^{2}}}},{a_{0402} = {\frac{K_{p\; 2}}{T_{v}} - T_{i\; 2}}},{a_{0404} = {{{- \frac{K_{p\; 2}\gamma\; u_{s\; 0}L_{m}}{\alpha\; T_{j}S_{B}L_{s}}}a_{0505}} = {- \frac{1}{T_{x}}}}},{a_{0511} = \frac{1}{T_{x}}},{a_{0605} = {\frac{K_{pcvc}}{T_{x}} - T_{icvc}}},{a_{0611} = {- \frac{K_{pcvc}}{T_{x}}}},{a_{0705} = {K_{pssc}\left( {\frac{K_{pcvc}}{T_{x}} - T_{icvc}} \right)}},{a_{0706} = T_{issc}},{a_{0707} = {- \frac{K_{pssc}}{L}}},{a_{0708} = {\frac{K_{pssc}R}{L} - T_{issc}}},{a_{0711} = {- \frac{K_{pssc}K_{pcvc}}{T_{x}}}},{a_{0807} = \frac{1}{L}},{a_{0808} = {- \frac{R}{L}}},{a_{0909} = {- \frac{K_{pssc}}{L}}},{a_{0910} = {\frac{K_{pssc}R}{L} - T_{issc}}},{a_{1009} = \frac{1}{L}},{a_{1010} = {- \frac{R}{L}}},{a_{1101} = \frac{\gamma\; u_{s\; 0}L_{m}i_{{qr}\; 0}}{\alpha\;{Cu}_{{dc}\; 0}L_{s}}},{a_{1104} = \frac{\gamma\; u_{s\; 0}{L_{m}\left( {1 - \omega_{m\; 0}} \right)}}{\alpha\;{Cu}_{{dc}\; 0}L_{s}}},{a_{1108} = \frac{\gamma\; u_{s\; 0}}{\;{Cu}_{{dc}\; 0}}},{a_{1111} = {- \frac{\gamma\;{u_{s\; 0}\left\lbrack {i_{1d\; 0} - \frac{\left( {1 - \omega_{m\; 0}} \right)L_{m}i_{{qr}\; 0}}{\alpha\; L_{s}}} \right\rbrack}}{{Cu}_{{dc}\; 0}^{2}}}}$ $B = \overset{\begin{matrix} {{\mspace{101mu}\mspace{31mu}}{\Delta\; u_{s}}} & {\Delta\theta} & \; \end{matrix}}{\begin{bmatrix} {\Delta{\overset{.}{\omega}}_{m}^{*}} & b_{0101} & 0 \\ {\Delta{\overset{.}{\omega}}_{ref}^{*}} & 0 & 0 \\ {\Delta{\overset{.}{i}}_{dr}} & b_{0301} & 0 \\ {\Delta{\overset{.}{i}}_{qr}} & b_{0401} & 0 \\ {\Delta{\overset{.}{x}}_{1}} & 0 & 0 \\ {\Delta{\overset{.}{i}}_{1{dref}}} & 0 & 0 \\ {\Delta{\overset{.}{u}}_{1{dref}}} & 0 & 0 \\ {\Delta{\overset{.}{i}}_{1d}} & 0 & 0 \\ {\Delta{\overset{.}{u}}_{1{dref}}} & 0 & 0 \\ {\Delta{\overset{.}{i}}_{1q}} & 0 & 0 \\ {\Delta{\overset{.}{u}}_{dc}} & b_{1101} & 0 \end{bmatrix}}$ ${b_{0101} = {- \frac{\gamma\; L_{m}i_{{qr}\; 0}}{\alpha\; L_{s}S_{B}}}},{b_{0401} = \frac{\gamma\; K_{p\; 2}L_{m}i_{{qr}\; 0}}{\alpha\; T_{j}L_{s}S_{B}}},{b_{1101} = \frac{\gamma\left\lbrack {i_{1\; d\; 0} - \frac{\left( {1 - \omega_{m\; 0}} \right)L_{m}i_{{qr}\; 0}}{\alpha\; L_{s}}} \right\rbrack}{{Cu}_{{dc}\; 0}}},{b_{0301} = {{\frac{L_{s}T_{i\; 1}K_{p\; 1}\gamma\; L_{m}\alpha}{\left( {{\alpha\; L_{s}} + {K_{p\; 1}\gamma\; u_{s\; 0}L_{m}}} \right)^{2}}\left\lbrack {Q_{ref} - {\frac{\gamma\; u_{s\; 0}}{L_{s}}\left( {\frac{i_{{dr}\; 0}L_{m}}{\alpha} - \frac{u_{s\; 0}}{\omega_{s}}} \right)}} \right\rbrack} + {\frac{\alpha\; T_{i\; 1}}{{\alpha\; L_{s}} + {K_{p\; 1}\gamma\; u_{s\; 0}L_{m}}}\left( {\frac{2\gamma\; u_{s\; 0}}{\omega_{s}} - \frac{\gamma\; i_{{dr}\; 0}L_{m}}{\alpha}} \right)}}}$ $\begin{matrix} \; \\ {C =} \end{matrix}\begin{matrix} \begin{matrix} {\mspace{95mu}{\Delta\;\omega_{m}^{\prime}}} & {\Delta\;\omega_{ref}^{\prime}} & {\Delta\; i_{dr}\begin{matrix} {\mspace{20mu}{\Delta\; i_{qr}}} & {\Delta\; x_{1}} & {\Delta\; i_{1{dref}}} & {\Delta\; u_{1{dref}}} & {\Delta\; i_{1d}} & {\Delta\; u_{1{dref}}} & {\Delta\; i_{1q}} \end{matrix}} & {\Delta\; u_{d}} \end{matrix} \\ {\begin{matrix} {\Delta\; P_{g}} \\ {\Delta\; Q_{g}} \end{matrix}\begin{bmatrix} 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & c_{0104} & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & c_{0108} & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 \\ 0 & \mspace{11mu} & 0 & \mspace{11mu} & c_{0203} & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & 0 & \mspace{11mu} & c_{0210} & \mspace{11mu} & 0 \end{bmatrix}} \end{matrix}$ ${c_{0104} = \frac{\gamma\; u_{s\; 0}L_{m}}{\alpha\; L_{s}}},{c_{0108} = {{- \gamma}\; u_{s\; 0}}},{c_{0203} = \frac{\gamma\; u_{s\; 0}L_{m}}{\alpha\; L_{m}}},{c_{0210} = {\gamma\; u_{s\; 0}}}$ $\begin{matrix} \; \\ {D =} \end{matrix}\begin{matrix} \begin{matrix} {\mspace{85mu}{\Delta\; u_{s}}} & {\Delta\;\theta} \end{matrix} \\ \begin{matrix} \begin{matrix} {\Delta\; P_{g}} \\ {\Delta\; Q_{g}} \end{matrix} & \begin{bmatrix} d_{0101} & 0 \\ d_{0201} & 0 \end{bmatrix} \end{matrix} \end{matrix}$ ${d_{0101} = {\gamma\left( {\frac{L_{m}i_{{qr}\; 0}}{\alpha\; L_{s}} - i_{1d\; 0}} \right)}},{d_{0201} = {\gamma\left( {i_{1\; q\; 0} - \frac{{2\alpha\; u_{s\; 0}} - {\omega_{s}L_{m}i_{{dr}\; 0}}}{\alpha\;\omega_{s}L_{s}}} \right)}}$ E = 0_(2 × 11) $\begin{matrix} \; \\ {F =} \end{matrix}\begin{matrix} \begin{matrix} {\mspace{85mu}{\Delta\; u_{s}}\;} & {\mspace{20mu}{\Delta\;\theta}\mspace{11mu}} \end{matrix} \\ \begin{matrix} \begin{matrix} {\Delta\; P_{g}} \\ {\Delta\; Q_{g}} \end{matrix} & \begin{bmatrix} f_{0101} & f_{0102} \\ f_{0201} & f_{0202} \end{bmatrix} \end{matrix} \end{matrix}$ f₀₁₀₁ = −m[E(G_(line)cos  θ₀ + B_(line)sin  θ₀) + 2 mu_(s 0)G_(line)] f₀₁₀₂ = mu_(s 0)E(G_(line)sin  θ₀ − B_(line)cos  θ₀) f₀₂₀₁ = −m[E(G_(line)sin  θ₀ − B_(line)cos  θ₀) + 2mu_(s 0)B_(line)] f₀₂₀₂ = −mu_(s 0)E(G_(line)cos  θ₀ + B_(line)sin  θ₀) $m = {\frac{k_{T}U}{u_{s\; 0}} = {k_{T}\sqrt{\frac{3}{2}}}}$ wherein, m is the ratio between effective value of high-side voltage of transformer and phase voltage amplitude of a generator stator. k_(T) is the k_(T) boosting transformer ratio.

A part of system parameters of this simulation model are shown in the following table.

TABLE 1 System Parameter Table Value/ Parameter Value/unit Parameter Value/unit Parameter unit R_(s) 0.0054 p.u. R_(r) 0.00607 p.u. L_(sσ) 0.1 p.u. L_(rσ) 0.11 p.u. L_(m) 4.5 p.u. ω_(s) 314.15 rad/s D 0.0001 p.u. A_(blade) 5026.55 m2 P 1.225 kg/m3 C_(p) 0.28 C 0.0078 F L 1 mH K_(pssc) 1 p.u. T_(issc) 0.1 p.u. T_(x) 0.04 μs

Analysis results of eigenvalues of system matrix A′ should refer to Table 2.

TABLE 2 Analysis Results of Eigenvalues of A′ Oscillation Most Second Damping frequency/ relevant relevant No. Eigenvalue ratio Hz variable variable λ₁ −1214.3316 1 0 u_(1dref) i_(1d) λ_(2,3) −4.8122 ± 130.8084i 0.0368 20.8188 u_(dc) i_(1d) λ₄ −1199.9167 1 0 u_(1qref) i_(1q) λ_(5,6) −0.7275 ± 1.0064i  0.5858 0.1602 ω*_(m) i_(qr) λ₇ −0.0010 1 0 x₁ λ₈ −0.1000 1 0 i_(1dref) λ₉ −0.0833 1 0 i_(1q) u_(1qref) λ₁₀ −20 1 0 ω*_(ref) λ₁₁ −0.8193 1 0 i_(dr)

From the analysis results of Table 2, all the eigenvalues of the system matrix A′ have negative real part, which represents the system is of small disturbance stability at this steady-state operation point. However, upon careful observation, it can figure out that the eigenvalues that correspond to evanescent modes λ₇, λ₈ and λ₉ are close to origin, so that unstability is easy to happen. Besides, damping of oscillation mode λ_(2,3) is small, which is only 0.0368, so that it is weak damping mode. For the above modes, it can figure out that small disturbance stability margin of this system is relatively low and when steady-state operating point shifts, unstability might happen.

After finding out the mode needs to be focused, taking Step 2 to solve sensitivities of the modes λ_(2,3), and λ₇, λ₈ and λ₉ to the nonzero elements in A′. For the expression of elements, the state matrix A can be used instead. Detailed results should refer to Table 3.

TABLE 3 Sensitivity of Weak Damping Mode to the Elements in A′ Corresponding Mode Element expression in A Sensitivity Abandon or not λ_(2,3) a₀₇₀₆′ T_(issc) 0.4950 − 0.0716i No a₀₇₀₈′ $\frac{K_{pssc}R}{L} - T_{issc}$ 0.4028 − 0.1041i No a₁₁₀₆′ 0 −0.0013 + 0.9443i  Yes λ₇ a₀₅₀₅′ $- \frac{1}{T_{x}}$  1.0011 No a₀₅₁₁′ $\frac{1}{T_{x}}$  1.0011 No λ₈ a₀₆₀₆′ 0  1.0013 Yes a₀₇₀₆′ T_(issc) −1.0015 No a₀₈₀₆′ 0 −1.0014 Yes λ₉ a₀₉₀₉′ $- \frac{K_{pssc}}{L}$  0.8334 No a₀₉₁₀′ $\frac{K_{pssc}R}{L} - T_{issc}$ −0.8334 No

It can be concluded from Table 3 that T_(issc), K_(pssc), T_(x), L and R are relatively more sensitive to the weak damping mode happened to this system. Conditions about the changes of different system modes are researched below when the above parameters change (the initial parameter values should refer to Table 1):

1) Conditions about the Changes of System Eigenvalues when T_(issc) Changes

Taking the change interval of T_(issc) as [0, 10], selecting parameter nodes as T_(issc)=0.1, 1, 5, 10 and observing simulation steady-state solution. It is concluded from the results that when T_(issc) changes, simulation steady-state solution is basically unchanged. Take ΔT_(issc)=1, change track of system eigenvalues should refer to FIGS. 2a and 2 b.

It can be concluded from FIGS. 2a and 2b that when T_(issc)=0, the system will have a zero eigenvalue. With T_(issc) increases, damping ratio λ_(2,3) of oscillation mode will decrease slightly and damping of evanescent modes λ₈ and λ₉ will increase obviously. As when T_(issc) increases, damping ratio of oscillation mode λ_(2,3) decreases little, but damping of evanescent modes λ₈ and λ₉ increases obviously, so that T_(issc) can be increased properly to improve small disturbance stability of the system, taking {circumflex over (T)}_(issc)=10.

2) Conditions about the Changes of System Eigenvalues when K_(pssc) Changes

Replacing original T_(issc) with {circumflex over (T)}_(issc)=10 obtained from 1), setting the change interval of K_(pssc) as [0, 10], selecting parameter nodes as K_(pssc)=0.1, 1, 5, 10 and observing simulation steady-state solution. It is concluded from the results that when K_(pssc) changes, simulation steady-state solution is basically unchanged. Setting ΔK_(pssc)−1, the changing track of the system eigenvalues should refer to FIGS. 3a and 3 b.

It can be concluded from FIGS. 3a and 3b that when K_(pssc)=0, positive real part happens to oscillation mode λ_(2,3)(λ_(2,3)=35.0110±16.4581i), the steady-state of the system becomes unstable and evanescent modes λ₁ and λ₄ become a group of oscillation modes. With K_(pssc) increases, damping of oscillation mode λ_(2,3) will increase, but the eigenvalues that correspond to weak evanescent modes λ₈ and λ₉ move from negative real axis to origin point and its degree of stability decreases. Thus, if K_(pssc) is too big or too small, the degree of stability of the system will be decreased. It can be figured out from FIGS. 3a and 3b that when K_(pssc)>3, damping of mode λ_(2,3) will increase slowly, but λ₈ and λ₉ still move to the origin point, so that take {circumflex over (K)}_(pssc)=2.

3) Conditions about the Changes of System Eigenvalues when T_(x) Changes

It can be concluded from Table 3 that T_(x) is only highly sensitive to λ₇. T_(x) is an inertia time constant while measuring DC voltage, the value of which is small. Setting the change interval as [0.01,0.1], it can be concluded from the results that when T_(x) changes, the simulation stability solution is basically unchanged. When ΔT_(x)=0.01, from simulation results, it figures out that the system eigenvalue that changes T_(x) is also basically unchanged, which indicates that weak damping of mode λ₇ is unrelated to system parameters and it is caused by system structure. Thus, if value of T_(x) does not change, keeping T_(x)=0.04.

4) Conditions about the Changes of System Eigenvalues when L Changes

L is the equivalent conversion inductance of a grid-side converter, setting change interval of L as [0.5,2] mH, selecting parameter nodes as L=0.5,1,1.5,2 and observing simulation steady-state solution. It is concluded from the results that when K_(pssc) changes, simulation steady-state solution is basically unchanged. Setting ΔL=0.1 mH, the change track of system eigenvalues should refer to FIGS. 4a and 4 b.

It can be concluded from FIGS. 4a and 4b that with L increases, damping of oscillation mode λ_(2,3) will decrease obviously and damping ratio also reduces; eigenvalue that corresponds to evanescent mode λ₉ tends to move away from an imaginary axis, but it is not obvious and can be ignored; strong evanescent modes λ₁ and λ₄ move towards the origin rapidly and damping decreases, but it still does not belong to strong damping and influence little on system stability. Therefore, if it is allowed, L is smaller, the stability margin of the system is higher. Take {circumflex over (L)}=0.5 mH here.

5) Conditions about the Changes of System Eigenvalues when R Changes

R is the AC side equivalent resistance of a grid-side converter. Upon analysis, it indicates that when R>0.5Ω, oscillation will happen to the system; when R<0.5Ω, value of R has very little influence on system eigenvalues. On the basis of the above conclusions, take {circumflex over (R)}=0.2Ω.

Until now, optimization for the parameters that correspond to the four concerned modes is completed. Table 4 is about the comparison between key mode eigenvalues of initial parameters and optimal parameters of the system. It can be concluded through analyzing the above contents, T_(x) and R influence little on eigenvalues. Thus, in Table 4, initial parameters and optimal parameters are set as T_(x)=0.04 and R=0.2Ω.

TABLE 4 Comparison between Key Mode Eigenvalues of Initial Parameters and Optimal Parameters Parameter Initial parameter Optimal parameter Mode T_(issc) = 0.1 {circumflex over (K)}_(pssc) = 1 {circumflex over (L)} = 1 mH {circumflex over (T)}_(issc) = 10 {circumflex over (K)}_(pssc) = 2 {circumflex over (L)} = 0.5 mH λ_(2,3) −4.8122 ± 130.8084i −9.6379 ± 137.1746i λ₇ −0.0010 −0.0010 λ₈ −0.1000 −5.0028 λ₉ −0.0833 −4.5502

It can be concluded from Table 4 that except mode λ₇, damping of other dominant modes has been raised for a certain extent. For λ₇, its weak damping is not caused for selecting improper parameters, so that it needs to be improved through changing system control structure or adding other control devices; for oscillation λ_(2,3), its damping has been raised from 0.0368 to 0.0701; λ₈ and λ₉ have changed from weak damping mode to strong damping mode and small disturbance stability margin of the system has been improved obviously.

What is said above has totally verified that while researching the weak damping mode generated after double-fed unit gets access to power grid, system damping can be effectively improved without adding other control means through optimizing system parameters with sensitivity analysis. Compared with traditional optimization methods. The present invention will greatly decrease unnecessary calculated amount and improve efficiency.

This test system is only a preferred embodiment of the present invention, but protection scope of the present invention is not confined to the embodiment. Any changes or replacements that the person of ordinary skill in the art can figure out easily within the technical range disclosed by the present invention should be included in the protection scope of the present invention. 

The invention claimed is:
 1. A computer implemented method for improving small disturbance stability after a double-fed power generation unit gets access to a power transmission system characterized in that the method comprises the following steps executed on a processor: step 1: building complete mathematical models for the double-fed power generation unit, the mathematical models mainly include an aerodynamic model, a power generator model, a mechanical model and a control system model; listing a system state equation and an output equation and then building a small disturbance mathematical model Δ{dot over (x)}=A′Δx by integrating power flow equation of the system after double-fed power generation unit gets access to the system; step 2: calculating a left modal matrix ψ and a right modal matrix φ of a matrix A′, determining the sensitivity of unstable modes or weak damping modes to matrix A′ with the formula $\frac{\partial\lambda_{i}}{\partial a_{k\; j}} = {\psi_{ik}\phi_{ji}}$  and finding out one to two nonzero elements a′_(ij) with the highest sensitivity in the matrix; analysis indicating that at a low and middle frequency band concerned by small disturbance stability of the system, the difference between eigenvalue of the matrix A′ and that of the state matrix A is not very large and as expression of A′ is very complex, element a_(ij) in A is used to make sensitivity analysis instead of the element a′_(ij) in A′; step 3: for controller or system parameters that can be set in a_(ij), changing value of these parameters within a certain interval, observing steady-state value of the variable required while calculating eigenvalue of matrix A′ in simulation results, then putting the steady-state value in A′ to solve the eigenvalues that correspond to each group of parameters and drawing a chart on the change track of eigenvalues; if eigenvalues are disperse, a part of the overlapped eigenvalues need to be locally enlarged to observe the trend for eigenvalues of dominant modes of the system; step 4: if there are other parameters that can be set in a_(ij), then repeating Step 3; step 5: comprehensively analyzing the chart on modal eigenvalues change with the track of parameters in Step 4, adjusting the parameters in Step 4, then selecting appropriate parameter combination upon comparison, with which both dominant modal damping and small disturbance stability margin of the system can be obviously improved.
 2. The computer implemented method for improving small disturbance stability after double-fed power generation unit gets access to the power transmission system according to claim 1 characterized in that the system matrix A′ is built in the following methods: selecting an appropriate state variable, an input variable and an output variable, state equation, output equation and power flow equation of the system can be expressed in the following general forms: {dot over (x)}=f(x,u) y=g(x,u)   (1), y =h(x,u) wherein x is a state variable matrix, u is an input variable matrix, y is an output matrix; wherein linearizing an equation (1) at a steady-state operating point, it can be concluded as: Δ{dot over (x)}=AΔx+BΔu Δy=CΔx+DΔu   (2), Δy=EΔx+FΔu wherein, $\begin{matrix} {{A = {{\begin{bmatrix} \frac{\partial f_{1}}{\partial x_{1}} & \ldots & \frac{\partial f_{1}}{\partial x_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial f_{n}}{\partial x_{1}} & \ldots & \frac{\partial f_{n}}{\partial x_{n}} \end{bmatrix} B} = {{\begin{bmatrix} \frac{\partial f_{1}}{\partial u_{1}} & \ldots & \frac{\partial f_{1}}{\partial u_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial f_{n}}{\partial u_{1}} & \ldots & \frac{\partial f_{n}}{\partial u_{n}} \end{bmatrix} C} = \left\lbrack \begin{matrix} \frac{\partial g_{1}}{\partial x_{1}} & \ldots & \frac{\partial g_{1}}{\partial x_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial g_{n}}{\partial x_{1}} & \ldots & \frac{\partial g_{n}}{\partial x_{n}} \end{matrix} \right\rbrack}}}{D = {{\left\lbrack \begin{matrix} \frac{\partial g_{1}}{\partial u_{1}} & \ldots & \frac{\partial g_{1}}{\partial u_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial g_{n}}{\partial u_{1}} & \ldots & \frac{\partial g_{n}}{\partial u_{n}} \end{matrix} \right\rbrack E} = {{\left\lbrack \begin{matrix} \frac{\partial h_{1}}{\partial x_{1}} & \ldots & \frac{\partial h_{1}}{\partial x_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial h_{n}}{\partial x_{1}} & \ldots & \frac{\partial h_{n}}{\partial x_{n}} \end{matrix} \right\rbrack F} = \left\lbrack \begin{matrix} \frac{\partial h_{1}}{\partial u_{1}} & \ldots & \frac{\partial h_{1}}{\partial u_{n}} \\ \ldots & \ldots & \ldots \\ \frac{\partial h_{n}}{\partial u_{1}} & \ldots & \frac{\partial h_{n}}{\partial u_{n}} \end{matrix} \right\rbrack}}}} & (20) \end{matrix}$ joining with the Equation (2), it can be concluded as: Δ{dot over (x)}=A′Δx  (4) wherein, A′=A+B(F−D)⁻¹(C−E)  (5).
 3. The computer implemented method for improving small disturbance stability after double-fed power generation unit gets access to the power transmission system comprising claim 1 characterized in that, in the Step 2, on the basis of the system matrix A′ obtained from Step 1, finding out the mode λ_(i),i=1,2, . . .,m that needs to be focused, wherein “m” refers to the number of unstable or weak damping modes; then finding out left modal matrix ψ′ and right modal matrix φ′ of matrix A′: for the left modal matrix ψ′: ψ′=[ψ′₁ ^(T)ψ′₂ ^(T) . . . ψ′_(n) ^(T)]^(T)  (6), wherein, ψ′_(i)A′=λ_(i)ψ′_(i),i=1,2, . . . ,n  (7) “n” is the number of state variables; for right modal matrix φ′: φ′=[φ′₁φ′₂ . . . φ′_(n)]  (8) wherein, A′φ′_(i)=λ_(i)φ′_(i),i=1,2, . . . ,n  (9) sensitivity of eigenvalue λ_(i) to the element of A′ can be expressed as: $\begin{matrix} {\frac{\partial\lambda_{i}}{\partial a_{kj}^{\prime}} = {{\psi_{i}^{\prime}\frac{\partial A^{\prime}}{\partial a_{kj}^{\prime}}\phi_{i}^{\prime}} = {\psi_{ik}^{\prime}\phi_{ji}^{\prime}}}} & (10) \end{matrix}$ the sensitivity of eigenvalue λ_(i) to a′_(kj) quantizes the change scope of λ_(i) when a′_(kj) changes, namely, when a′_(kj) changes, $\frac{\partial\lambda_{i}}{\partial a_{kj}^{\prime}}$  is larger, λ_(i) changes more obviously; thus after getting eigenvalues of A′, for the unstable mode, weak damping evanescent mode and weak damping ratio oscillation mode that may influence small disturbance stability of the system directly, the nonzero element that is most sensitive to this mode can be found according to the above method; adapting an Equation (5) into: A′=A+A _(other)  (11) wherein, A _(other) =B(F−D)⁻¹(C−E)  (12) from Equation (11), it can be visually seen that state matrix A is a component of system matrix A′, so that corresponding elements of A also exist in A′; then after obtaining $\frac{\partial\lambda_{i}}{\partial a_{kj}^{\prime}},$  i=1,2, . . . ,m and finding out set of the elements that are highly sensitive to mode λ_(i),i=1,2, . . . ,m, analyzing with the component {a_(kj)} of state matrix that shares the same code with the elements in {a′_(kj)}; the reason for finding nonzero element is that: for a system with fixed structure, structure of its system matrix A′ is also fixed; if a_(kj)=0, no matter how to change parameters, a_(kj) remains unchanged.
 4. The computer implemented method for improving small disturbance stability after double-fed power generation unit gets access to the power transmission system according to claim 1 characterized in that, in the Step 3, on the basis of finding out the set {a_(kj)} of elements with high sensitivity in Step 2, finding out the adjustable controller parameters or system parameters in {a_(kj)}, wherein, k_(i),i=1,2, . . . ,t, “t” refers to the number of adjustable variables in {a_(kj)}.
 5. The computer implemented method for improving small disturbance stability after double-fed power generation unit gets access to the power transmission system according to claim 4 characterized in that, the k₁ refers that: letting k₁ change within set interval [a,b] and selecting several parameter nodes within this interval, then cycling calculate eigenvalues of A′ while k₁ is changing and selecting the optimal one {circumflex over (k)}₁.
 6. The computer implemented method for improving small disturbance stability after double-fed power generation unit gets access to the power transmission system according to claim 1 characterized in that, in the Step 4, on the basis of selecting {circumflex over (k)}₁, repeating Step 3 for other adjustable parameters k_(i),i=2, . . . ,t, until all the adjustable parameters are set.
 7. The computer implemented method for improving small disturbance stability after double-fed power generation unit gets access to the power transmission system according to claim 1 characterized in that, in the Step 5, comprehensively comparing analysis results of the eigenvalues of optimal parameter set {circumflex over (k)}_(i),i=1,2,, . . . ,t and original parameter set k_(i),i=1,2, . . . ,t; upon analysis, modal damping of the system can be greatly improved and small disturbance stability margin of the system can be intensified after optimizing parameters with sensitivity analysis. 